    
    module routines_specific
    
      use globalvalues
      use precisions
      use routines_generic
      use types_generic

      implicit none

    contains
    
    subroutine readInParms(readParams, numParamsComb,allInputsAllPeeps)
            implicit none
            
            integer, intent(in) :: numParamsComb
            type(typeReadInParams), intent(out) :: readParams(numParamsComb)
            type(typeAllInputs), intent(out) :: allInputsAllPeeps(numPeeps)
            character (130) :: path_in
            character (130) :: path_out
            character (130) :: in_params, in_inputs, in_inputsUnscaled
            integer ::  ixParamComb, ixElsaagent, ix
            real (kind=rk) :: readindata(8) ! This is number of parameters
            real (kind=rk) :: allInputs(61), allInputsUnscaled(61)
         
            
            call getpaths(path_in)
            
           in_params = trim(path_in) // 'inputParamsMain.txt'

            
            
            open (unit=110, file=in_params, action='read')
  
            !For programme
            do ixParamComb = 1, numParamsComb
                read (110, *) readindata
                readParams(ixParamComb)%beta = readindata(1)              
                readParams(ixParamComb)%gamma = readindata(2)
                readParams(ixParamComb)%obj = readinlogical(readindata(3), 3)     
                readParams(ixParamComb)%addMinLoad = readinlogical(readindata(4), 4)        
                readParams(ixParamComb)%housing = readinlogical(readindata(5), 5)     
                readParams(ixParamComb)%annOn  = readinlogical(readindata(6), 6)  
                readParams(ixParamComb)%unscaled     = readinlogical(readindata(7), 7)  
                readParams(ixParamComb)%discreteAnn   = readinlogical(readindata(8), 8)    
            end do
            

                in_inputs = trim(path_in) // 'annInputs.txt'
                in_inputsUnscaled = trim(path_in) // 'annInputsUnscaled.txt'        
            
            open (unit=111, file=in_inputs, action='read')
            open (unit=112, file=in_inputsUnscaled, action='read')
            
            do ixElsaagent = 1, numPeeps
                read(111,*) allInputs
                read(112,*) allInputsUnscaled
                allInputsAllPeeps(ixElsaagent)%weibullK = allInputs(1)
                allInputsAllPeeps(ixElsaagent)%weibullLam = allInputs(2)

                
            allInputsAllPeeps(ixElsaagent)%age = allInputs(3)
            allInputsAllPeeps(ixElsaagent)%objectiveAnn = allInputs(4)
            allInputsAllPeeps(ixElsaagent)%id = allInputs(5)
            allInputsAllPeeps(ixElsaagent)%statePen = allInputs(6)
            allInputsAllPeeps(ixElsaagent)%Assets = allInputs(7)
            
                do ix = 1, 51
                    allInputsAllPeeps(ixElsaagent)%hazards(ix) = allInputs(8+ ix)         
                    allInputsAllPeeps(ixElsaagent)%hazardsUnscaled(ix) = allInputsUnscaled(8+ ix)         
                end do
    
                allInputsAllPeeps(ixElsaagent)%imputedRent = allInputs(61)
            end do
            
            
            
    end subroutine readInParms
    !----------------------------------------------------------------!
    !----------------------------------------------------------------!
        
    real (kind=rk) function felicity(c, params)
        implicit none

! Arguments
        real (kind=rk), intent (in) :: c
        type (typeParams), intent (in) :: params
        
        real (kind=rk) :: consWithFree
      
        
        consWithFree =  (c + params%freeMinCons) + params%imputedRentForUtil 
        ! Check
        if (consWithFree.lt.0.0_rk) then
            write(*,*) 'error in felicity function, consumption is negative'
            write(*,*) c
            stop
        end if
                
            if ((params%gamma.gt.(1.0_rk-1.0d-2)).and.(params%gamma.lt.(1.0_rk+1.0d-2))) then
              felicity =    log(consWithFree)
            else
              felicity =   ( ((consWithFree**params%OneMinusGamma)/params%OneMinusGamma) )
            end if
            


            
    end function felicity
    
    

    
    
    subroutine getgrids(tech,grids)
        implicit none
        type(typeTech),  intent(in) :: tech  
        type(typegrids), intent(inout) :: grids
  
        real (kind=rk), dimension(grids%maxT+1) :: allTAssMins, allTAssMax
        real (kind=rk), dimension(grids%maxT+1) :: allTYMins, allTYMax, allTYBreak1, allTYBreak2
        integer :: ixT, ixY
        

        do ixT = 1, grids%maxT+1
            allTAssMins(ixT) = 0.0_rk !grids%assMin  
            allTAssMax(ixT) =  max(grids%assMax, grids%assMax/1.05_rk +1.0_rk)  
        end do        

        

     do ixT = 1,  grids%maxT+1 
     call getnodes(grids%A(:, ixT), 0.0_rk, grids%assMax/1.05_rk , grids%maxA, 0.05_rk) 
     end do    
     
      
    
     !! Get consumption grid for optimization
        call getnodes(grids%c, 0.0_rk, 0.99999_rk, grids%cPoints, 0.00_rk) 
       
     !! Get annuitization grid for optimaztion
        call getnodes(grids%ann, 0.0_rk, 0.99999_rk, grids%annPoints, 0.00_rk) 
       
       do ixT = tech%retireAge,  grids%maxT+1 
            do ixY = 1, grids%maxYO,1
                grids%YOld(ixY, ixT) = (grids%assMax/1.05_rk) * tech%annRates(1)  *grids%ann(ixY)
            end do
        end do
 
                
    end subroutine getgrids 
    

    
    
    subroutine printTLStuff(TLExpDetails,singleRun, paramComb, params)
              
        use globalvalues !This is ONLY to allow getpaths to run

        implicit none
        real (kind=rk) :: TLExpDetails(4, 4, numpeeps) 
        logical, intent(in) :: singleRun
        integer, intent(in) :: paramComb
        type (typeParams), intent (in) :: params  
        
        character (10) :: tempformat_numpeepcols
        character (20) :: format_numpeepcols
        character (20) :: format_numpeepcols_int
        character (130) :: path_in
        character (130) :: path_out
        character (130) :: outfile

        real (kind=rk) :: TLval1(4, numpeeps)
        real (kind=rk) :: TLval2(4, numpeeps)
        real (kind=rk) :: TLann(4, numpeeps)
        real (kind=rk) :: TLass(4, numpeeps)
        real (kind=rk) :: mult
        
        integer :: requiredl, ios
       
     
! Opening files for output
        call getpaths(path_in, path_out,  paramComb)
            
        if ((params%gamma).gt.1) then
            mult = 1.0_rk/((10000.0_rk)**(1-params%gamma))
        else
            mult = 1.0_rk
        end if
        
        TLval1 = TLExpDetails(:,1,:) * mult
        TLval2 = TLExpDetails(:,2,:)  * mult
        TLann = TLExpDetails(:,3,:)
        TLass = TLExpDetails(:,4,:)
        
!  In printing it I have to tell format how many columns there are.
! If there are 10 columns and the format of each entry is f15.15, the format I give write needs to be 10f15.5.
! In the first row I want the to print time period 1 for each individual. So I want to read across the first
! row of sdfsim (numyears, howmanysims). But write(*,*) reads down each columns first. So I need to take
! the tranpose
! so the first column of the matrix is in fact the first period for  rather than all periods for person 1

! Finally, I need to generate a string that is 'x'f15.5 where 'x' is the number of periods. I do this below
        write (tempformat_numpeepcols, '(I10)') numPeeps
        format_numpeepcols = '(' // trim(adjustl(tempformat_numpeepcols)) // 'f20.8' // ')'
        format_numpeepcols_int = '(' // trim(adjustl(tempformat_numpeepcols)) // 'I5' // ')'
        
        outfile = trim(path_out) // 'TLval1.out'               
        !inquire (iolength=requiredl) transpose(TLval1)
        !open (unit=250, file=outfile, recl=requiredl, action='write')
        open (unit=250, file=outfile,action='write')
        !write (*, trim(format_numpeepcols),  IOSTAT=ios) transpose(TLval1)
        write (250, trim(format_numpeepcols),  IOSTAT=ios) transpose(TLval1)

        outfile = trim(path_out) // 'TLval2.out'               
        open (unit=250, file=outfile,action='write')
        write (250, trim(format_numpeepcols),  IOSTAT=ios) transpose(TLval2)
        
        outfile = trim(path_out) // 'TLann.out'               
        open (unit=250, file=outfile,action='write')
        write (250, trim(format_numpeepcols),  IOSTAT=ios) transpose(TLann)
        
        outfile = trim(path_out) // 'TLass.out'               
        open (unit=250, file=outfile,action='write')
        write (250, trim(format_numpeepcols),  IOSTAT=ios) transpose(TLass)      
        

        close(250)
        
  end subroutine printTLStuff
  
    
          subroutine printsimdetail(profilesAnn, sim,  paramComb)
              
        use globalvalues !This is ONLY to allow getpaths to run

        implicit none
        type(typeProfilesAnn), intent(in) :: profilesAnn
        type(typeSim), intent(in)  :: sim    
        integer, intent(in) :: paramComb

          
        character (10) :: tempformat_numpeepcols
        character (20) :: format_numpeepcols
        character (20) :: format_numpeepcols_int
        character (130) :: path_in
        character (130) :: path_out
        character (130) :: outfile
               
        integer :: requiredl, ios
       
     
! Opening files for output
        call getpaths(path_in, path_out,paramComb)
            
  if (printMicroData) then
 !In printing it I have to tell format how many columns there are.
! If there are 10 columns and the format of each entry is f15.15, the format I give write needs to be 10f15.5.
! In the first row I want the to print time period 1 for each individual. So I want to read across the first
! row of sdfsim (numyears, howmanysims). But write(*,*) reads down each columns first. So I need to take
! the tranpose
! so the first column of the matrix is in fact the first period for  rather than all periods for person 1

! Finally, I need to generate a string that is 'x'f15.5 where 'x' is the number of periods. I do this below
        write (tempformat_numpeepcols, '(I10)') numPeeps
        format_numpeepcols = '(' // trim(adjustl(tempformat_numpeepcols)) // 'f20.8' // ')'
        format_numpeepcols_int = '(' // trim(adjustl(tempformat_numpeepcols)) // 'I5' // ')'
        
        if (fullPrinting) then
            outfile = trim(path_out) // 'cons.out'               
            inquire (iolength=requiredl) transpose(profilesAnn%cons)
            open (unit=250, file=outfile, recl=requiredl, action='write')
            write (250, trim(format_numpeepcols),  IOSTAT=ios) transpose(profilesAnn%cons)
            close(250)
                                    
            outfile = trim(path_out) // 'assStart.out'               
            inquire (iolength=requiredl) transpose(profilesAnn%assStart)
            open (unit=250, file=outfile, recl=requiredl, action='write')
            write (250, trim(format_numpeepcols),  IOSTAT=ios) transpose(profilesAnn%assStart)
            close(250)
        
            outfile = trim(path_out) // 'assEnd.out'               
            inquire (iolength=requiredl) transpose(profilesAnn%assEnd)
            open (unit=250, file=outfile, recl=requiredl, action='write')
            write (250, trim(format_numpeepcols),  IOSTAT=ios) transpose(profilesAnn%assEnd)
            close(250)
        
            outfile = trim(path_out) // 'Y.out'               
            inquire (iolength=requiredl) transpose(profilesAnn%assEnd)
            open (unit=250, file=outfile, recl=requiredl, action='write')
            write (250, trim(format_numpeepcols),  IOSTAT=ios) transpose(profilesAnn%Y)
            close(250)
        end if


        outfile = trim(path_out) // 'annProp.out'               
        inquire (iolength=requiredl) transpose(profilesAnn%annProp)
        open (unit=250, file=outfile, recl=requiredl, action='write')
        write (250, trim(format_numpeepcols),  IOSTAT=ios) transpose(profilesAnn%annProp)
        close(250)
        


       outfile = trim(path_out) // 'val.out'               
       inquire (iolength=requiredl) transpose(profilesAnn%val)
       open (unit=250, file=outfile, recl=requiredl, action='write')
        write (250, trim(format_numpeepcols),  IOSTAT=ios) transpose(profilesAnn%val*1000.0_rk)
        close(250)                   
       
        
                      
  end if
  
  
        call dotheprinting(yearsGlob, path_out,'cons', profilesAnn%consMean) 
        call dotheprinting(yearsGlob, path_out,'assStart', profilesAnn%assStartMean) 
        call dotheprinting(yearsGlob, path_out,'assEnd', profilesAnn%assEndMean)
        call dotheprinting(yearsGlob, path_out,'Y', profilesAnn%YMean)  
        call dotheprinting(yearsGlob, path_out,'annProp', profilesAnn%annPropMean)

        
          end subroutine printsimdetail

          

          
    
        subroutine dotheprinting(numyears, path_out,filestub, forprinting_mean)
          implicit none 
          integer, intent(in) :: numyears
          real (kind=highprecintkind), intent (in), dimension (numyears) :: forprinting_mean
          character (130), intent (in) :: path_out
          character (*), intent (in) :: filestub
          character (130) :: dumout_mean, dumout_sdev
          integer :: requiredl

          dumout_mean = trim(path_out)  // trim(filestub) // 'mean.out'
          dumout_sdev = trim(path_out) // trim(filestub) // 'sdev.out'

          
          inquire (iolength=requiredl) forprinting_mean
          open (unit=300, file=dumout_mean, recl=requiredl, action='write')
          !write (300, trim(format_numtypecols)) forprinting_mean
          write (300,'(f15.5)') forprinting_mean
          close (300)

          !inquire (iolength=requiredl) forprinting_sdev
          !open (unit=300, file=dumout_sdev, recl=requiredl, action='write')
          !!write (300, trim(format_numtypecols)) forprinting_sdev
          !write (300, '(f15.5)') forprinting_sdev
          !close (300)


          !close (300)

        end subroutine dotheprinting
! -------------------------------------------
          
       !        subroutine concat(dumfilestub, dumout_mean, dumout_sdev, dumout_skew, dumout_kurt, dumout_min, dumout_max)
       !   implicit none

        !  character (100), intent (in) :: dumfilestub
        !  character (100), intent (out) :: dumout_mean
        !  character (100), intent (out) :: dumout_sdev
        !  character (100), intent (out) :: dumout_skew
        !  character (100), intent (out) :: dumout_kurt
        !  character (100), intent (out) :: dumout_min
        !  character (100), intent (out) :: dumout_max

        !  dumout_mean = trim(path_out) // trim(dumfilestub) // 'mean.out'
        !  dumout_sdev = trim(path_out) // trim(dumfilestub) // 'sdev.out'


        !end subroutine concat   
          

subroutine getannuityratesSingle(years,  annAdminLoad, r, currentAge,survObj, annRates)
  
  implicit none

  !type(typeTech),  intent(in) :: tech  
  integer, intent(in) :: years
  real (kind=rk), intent(in) :: survObj(years)
  real (kind=rk), intent(in) :: annAdminLoad, r
  integer, intent(in) :: currentAge
  real (kind=rk), intent(out) :: annRates(years)

  
  
  !For program
  real (kind=rk) :: capitalf(years) !This and the following follow the notation in my notes
  integer :: ixtoday
  integer :: tau, ageCounter
  real (kind=rk) :: temp
  real (kind=rk) :: inc
  
  capitalf(years) = 0.0_rk
  
  
  !Elsewhere in the programme - I interpret survprob(ixtoday) as the probability as surviving
  !until ixtoday + 1. So, notwithstanding the fact that in my notes I have written this object
  !as s_{ixtoday + 1) I follow the convention here. 
  
  
  ! Note compared to some of my other code, here the administrative load is e.g. 0.1 rather than e.g. 1.1
      
  
      do ixtoday = years-1, 1, -1
         tau = years - ixtoday !years left
  
           capitalf(ixtoday)  =   (survObj(ixtoday) * (1.0_rk + capitalf(ixtoday+1)) ) / (1 + r)        
           annRates(ixtoday) = (1.0_rk/( ((1.0_rk)* (1.0_rk)) *  capitalf(ixtoday))) *(1 -  annAdminLoad)    
      end do
  
annRates(years) = 999.0_rk
  
    
end subroutine getannuityratesSingle          


! I want to get subjectively fair annuity rates for printing. Issues with very low survival probabilities. I needed to make some edits. But to avoid unintentded consequences for the battle hardened routine
! above, I did a copy and paste

subroutine getannuityratesSingleForPrinting(years,  annAdminLoad, r, currentAge,survObj, annRates)
  
  implicit none

  !type(typeTech),  intent(in) :: tech  
  integer, intent(in) :: years
  real (kind=rk), intent(in) :: survObj(years)
  real (kind=rk), intent(in) :: annAdminLoad, r
  integer, intent(in) :: currentAge
  real (kind=rk), intent(out) :: annRates(years)

  
  
  !For program
  real (kind=rk) :: capitalf(years) !This and the following follow the notation in my notes
  integer :: ixtoday
  integer :: tau, ageCounter
  real (kind=rk) :: temp
  real (kind=rk) :: inc
  
  capitalf(years) = 0.0_rk
  
  
  !Elsewhere in the programme - I interpret survprob(ixtoday) as the probability as surviving
  !until ixtoday + 1. So, notwithstanding the fact that in my notes I have written this object
  !as s_{ixtoday + 1) I follow the convention here. 
  
  
  ! Note compared to some of my other code, here the administrative load is e.g. 0.1 rather than e.g. 1.1
      
  
      do ixtoday = years-1, 1, -1
         tau = years - ixtoday !years left
  
           capitalf(ixtoday)  =   (survObj(ixtoday) * (1.0_rk + capitalf(ixtoday+1)) ) / (1 + r)        

           if (capitalf(ixtoday)  .lt. 1d-20) then
               capitalf(ixtoday)  = 1d-20
           end if
           

            annRates(ixtoday) = (1.0_rk/( ((1.0_rk)* (1.0_rk)) *  capitalf(ixtoday))) *(1 -  annAdminLoad)    

      end do
  
annRates(years) = 999.0_rk
  
    
end subroutine getannuityratesSingleForPrinting

subroutine getData(grids, params, tech)
        type(typeGrids), intent (inout) :: grids
        type (typeParams), intent (inout) :: params
        type(typeTech),  intent(inout) :: tech  

        
        ! For program
        integer :: ixY, ix, ixAR1, ixT, ixTrans, ixElsaagent
        character (130) :: in_survprobs, in_detIncome, path_in, path_out
        
        real (kind=rk) :: totalSurvProb(52)
        !real (kind=rk) :: annRatesFair(tech%maxT), annRatesLoaded(tech%maxT)
        logical :: firstTime
        integer :: forPower
        
        call getpaths(path_in)                              
        
do ixElsaagent = 1, numpeeps        

        !! Survival Probabilities
        if (tech%survProbOn) then
  
          tech%survObj(1) = 1.0_rk
          tech%survObj(2:51) = 1-tech%allIn(ixElsaagent)%hazards(1:50)
          tech%survObj(52) = 0.0_rk
          
          tech%survObjUnscaled(1) = 1.0_rk
          tech%survObjUnscaled(2:51) = 1-tech%allIn(ixElsaagent)%hazardsUnscaled(1:50)
          tech%survObjUnscaled(52) = 0.0_rk
          
          ! For ages after 110 make surival probability tiny            
             do ix = 1,52   
              ! If age less than 109, check that survival probability is non-zero
                 !When age is 60 the first unit below gives 60 + 1 - 2 = 59. And the first object is the probability of surviving until 60 - the first period.                 
            if ((tech%allIn(ixElsaagent)%age + ix - 2).le.109) then !should have some probability of survival
                if ((tech%survObj(ix) .lt. 1d-20) .or. (tech%survObj(ix) .gt.1.0_rk)) then
                    write(*,*) 'error101' 
                    stop
                end if
            else if ((tech%allIn(ixElsaagent)%age + ix - 2).eq.110) then
                if ((tech%survObj(ix) .gt. 1d-20) .or. (tech%survObj(ix) .lt.0.0_rk)) then
                    write(*,*) 'error101' 
                    stop
                end if    
                tech%survObj(ix)  =  1d-20 ! if its zero it causes an error   
            else if ((tech%allIn(ixElsaagent)%age + ix - 2).gt.110) then
                if ((tech%survObj(ix) .gt. 1d-20) .or. (tech%survObj(ix) .lt.0.0_rk)) then
                    write(*,*) 'error101' 
                    stop
                end if    
                tech%survObj(ix)  =  1d-20 ! if its zero it causes an error                       
            end if

            if ((tech%allIn(ixElsaagent)%age + ix - 2).le.109) then !should have some probability of survival
                if ((tech%survObjUnscaled(ix) .lt. 1d-20) .or. (tech%survObjUnscaled(ix) .gt.1.0_rk)) then
                    write(*,*) 'error101' 
                    stop
                end if
            else if ((tech%allIn(ixElsaagent)%age + ix - 2).eq.110) then
                if ((tech%survObjUnscaled(ix) .gt. 1d-20) .or. (tech%survObjUnscaled(ix) .lt.0.0_rk)) then
                    write(*,*) 'error101' 
                    stop
                end if    
               tech%survObjUnscaled(ix)  =  1d-20 ! if its zero it causes an error   
            else if ((tech%allIn(ixElsaagent)%age + ix - 2).gt.110) then
                if ((tech%survObjUnscaled(ix) .gt. 1d-20) .or. (tech%survObjUnscaled(ix) .lt.0.0_rk)) then
                    write(*,*) 'error101' 
                    stop
                end if    
                tech%survObjUnscaled(ix)  =  1d-20 ! if its zero it causes an error                       
            end if
            
            end do

        else
            tech%survObj=1.0_rk    
            tech%survObjUnscaled = 1.0_rk
        end if
        
        firstTime = .true.
        do ix = 1,52   
             totalSurvProb(ix)  = exp( - ( ( real(ix-1, rk) /tech%allIn(ixElsaagent)%weibullLam)** tech%allIn(ixElsaagent)%weibullK)) 
             if (totalSurvProb(ix) .lt. 1d-10) then
                      totalSurvProb(ix)=  1d-10 ! / (10** (ix-forPower))
                 
             end if

 


             
             if (ix.eq.1) then
                 tech%surv(ix)  = totalSurvProb(ix)
             else if (totalSurvProb(ix-1).gt.1d-9) then
                 tech%surv(ix) = totalSurvProb(ix)/totalSurvProb(ix-1)
             else
                 tech%surv(ix) = 0.0_rk
             end if             
          
            
            if ((tech%allIn(ixElsaagent)%age + ix - 1).le.110) then !should have some probability of survival
             if (totalSurvProb(ix) .lt. 1d-20) then
                  tech%surv(ix)  =  1d-20 
             end if       
            else if ((tech%allIn(ixElsaagent)%age + ix - 1).gt.110) then
                tech%surv(ix)  =  1d-20 ! if its zero it causes an error
            end if
            
        end do
        
                  

        ! Get annuity rates
        if (tech%annOn) then
            call getannuityratesSingle(grids%maxT, 0.0_rk,  tech%r, tech%allIn(ixElsaagent)%age, tech%survObj, tech%annRatesAll(:,ixElsaagent, 1, 1)) !Fair Scaled
            call getannuityratesSingle(grids%maxT, tech%annAdminLoad, tech%r, tech%allIn(ixElsaagent)%age, tech%survObj, tech%annRatesAll(:,ixElsaagent, 2, 1)) !Loaded scaled
            call getannuityratesSingle(grids%maxT, tech%annAdminLoad, tech%r, tech%allIn(ixElsaagent)%age, tech%survObjUnscaled, tech%annRatesAll(:,ixElsaagent, 2, 2)) !Loaded unscaled
            call getannuityratesSingle(grids%maxT, 0.0_rk,  tech%r, tech%allIn(ixElsaagent)%age, tech%survObjUnscaled, tech%annRatesAll(:,ixElsaagent, 1, 2)) !Fair unscaled        
            
            ! I want to get subjectively fair annuity rates. Issues with very low survival probabilities. I needed to make some edits. But to avoid unintentded consequences for the battle hardened routine
            ! above, I did a copy and paste
            call getannuityratesSingleForPrinting(grids%maxT, 0.0_rk,  tech%r, tech%allIn(ixElsaagent)%age, tech%surv, tech%annRatesSubj(:,ixElsaagent)) !Subjective, just for Printing     
        else
            tech%annRates(:) = 0.0_rk
        end if 
        


        ! Now populate the things that change across people
        tech%survProbAll(:,ixElsaagent) = tech%surv
        tech%survProbAllObj(:,ixElsaagent) = tech%survObj            
        tech%survProbAllObjUnscaled(:,ixElsaagent) = tech%survObjUnscaled
end do



end subroutine getData

subroutine getMax(infoForMax, grids, valpol, params, tech, argmax, max)
        use globalvalues
        use types_generic
        implicit none


        
! Arguments
        type(typeGrids), intent(in) :: grids
        type(typeInfoForMax), intent(in) :: infoForMax
        type(typeValPol), intent(in) :: valpol
        type(typeParams), intent(in) :: params
        type(typeTech), intent(in) :: tech
        real (kind=rk), intent (out) :: argmax(6), max
        

! For program

        type(typeForOptRout) :: ForOptRoute 
        integer :: ixC, maxIx, maxIxCP(3), ixAnn 
        real (kind=rk) :: cons, gridOfCons(grids%cPoints), gridOfValues(grids%cPoints)
        real (kind=rk) :: gridOfAnnInc(grids%annPoints), gridOfValues_C(grids%cPoints)
        real (kind=rk) :: gridOfOptCons_byAPPC(grids%annPoints), gridOfValues_byAPPC(grids%annPoints) 
        
        integer :: maxIxAPPC(1) 
        real (kind=rk) :: AssetLeft, AnnPot 
        real (kind=rk) :: candidateArgMax(2), candidateMax

        real (kind=rk) :: resources 
            
! For the search for a local maximum
        real (kind=rk) :: lowerB, upperB


        
! In years in which there is an annuity choice
if (infoForMax%today.eq.(tech%retireAge - 1)) then


        !call setResourcesAndPCPaid(infoForMax, tech, 1, resources, ForOptRoute%PCpaidOpt) ! the 1 is ixPC
  resources = infoForMax%cashOnHand

      
          
                do ixAnn = 1, grids%annPoints 
                            if (((.not.(tech%annOn)) .and. (ixAnn .gt.1)).or. (tech%discreteAnn).and.(ixAnn.gt.1).and.(ixAnn.lt.grids%annPoints)) then
                                 gridOfValues_byAPPC(ixAnn) = -huge(rk)
                            else                        
                                    AnnPot  = grids%ann(ixAnn) * resources ! This is the pot of money that will be spend in an annuity
                                    gridOfAnnInc(ixAnn) = AnnPot * tech%annRates(infoForMax%today) ! This is the income flow that that pot will generate
                                    ForOptRoute%AssetForLiqOrCon = resources - AnnPot
                                    ForOptRoute%Y1 = 0 
                                    ForOptRoute%ixAnn = ixAnn  ! THIS IS USED TO SELECT OUT THE RELEVANT INCOME SLICE IN THE ONE DIMENSIONAL INTERPOLATION

                                                                          
                                                    do ixC = 1, grids%cPoints, 1       
                                                    cons =   ForOptRoute%AssetForLiqOrCon * grids%C(ixC)                              
                            
                                                    if (cons.lt.0.0_rk) then
                                                        write(*,*) 'Negative consumption in getMax'
                                                        stop
                                                    end if
                            
                            
                                                    gridOfCons(ixC) = cons
                                                    ForOptRoute%leftOverA =  (ForOptRoute%AssetForLiqOrCon - cons)*(1+tech%r)                                    
                                                    call getUtilPlusVal(infoForMax, ForOptRoute, grids, valpol, params, tech, cons, gridOfValues_C(ixC))                                 
                                                    end do 
                            
                                                maxIx = maxloc(gridOfValues_C, 1)     
                               
                                                ForOptRoute%annPot = AnnPot                    
                                                candidateMax = gridOfValues_C(maxIx)
                                                candidateArgMax(1) = gridOfCons(maxIx) 
                                                if ((doLocalAfterGrid).and.(infoForMax%today.ne.1)) then
                                                    call doLocal(gridOfOptCons_byAPPC(ixAnn), gridOfValues_byAPPC(ixAnn) ) 
                                                else
                                                    gridOfOptCons_byAPPC(ixAnn) = candidateArgMax(1)
                                                    gridOfValues_byAPPC(ixAnn) = candidateMax
                                                end if

                             end if                    
                    
                end do !ixAnn
                
  

            maxIxAPPC = maxloc(gridOfValues_byAPPC)                      
            max = gridOfValues_byAPPC(maxIxAPPC(1) )
            argmax(1) = gridOfOptCons_byAPPC(maxIxAPPC(1))
            argmax(2) = grids%ann(maxIxAPPC(1))

            

else    
        ForOptRoute%annPot = 0.0_rk ! The golden search will use this for assets next period, so needs to be zero here
    

        resources = infoForMax%cashOnHand
        
            if (resources.lt.0.0_rk) then ! can't afford to pay the participation charge
                write(*,*) 'resources less than zero in getMax'
                stop
            end if
            
                            
                             ForOptRoute%AssetForLiqOrCon = resources 
                                do ixC = 1, grids%cPoints, 1
                                       cons =  (ForOptRoute%AssetForLiqOrCon) * grids%C(ixC)  
                            
                                        if (cons.gt.ForOptRoute%AssetForLiqOrCon) then ! Shouldn't be possible but numerical issues mean that it is sometimes when consGrid == 1
                                            cons = ForOptRoute%AssetForLiqOrCon - 0.0000000001_rk
                                        end if               
                
                                        gridOfCons(ixC) = cons
                                        ForOptRoute%leftOverA =  (ForOptRoute%AssetForLiqOrCon - cons)*(1+tech%r)
                                        call getUtilPlusVal(infoForMax,  ForOptRoute, grids, valpol, params, tech, cons, gridOfValues_C(ixC))   
                                end do 
            
                               
                                maxIx = maxloc(gridOfValues_C, 1)     
                                if (doLocalAfterGrid) then
                                    candidateMax = gridOfValues_C(maxIx)
                                    candidateArgMax(1) = gridOfCons(maxIx) 
                                    call doLocal(argmax(1), max ) 
                                else
                                    argmax(1) = gridOfCons(maxIx) 
                                   max = gridOfValues_C(maxIx)
                                end if     
                    
            if (argmax(1) .lt. 0.0_rk) then
                write(*,*) 'negative consumption returned by getMax'
                stop
            end if
            


end if

contains


subroutine doLocal(argmaxDum, maxDum)
        implicit none
        real (kind=rk), intent(out) :: argmaxDum, maxDum
        
            ! Now implement the golden search
            if (maxIx.eq.1) then
                lowerB = gridOfCons(1)
                upperB = gridOfCons(2)
            else if (maxIx.eq.grids%cPoints) then 
                lowerB = gridOfCons(grids%cPoints-1)
                upperB = gridOfCons(grids%cPoints)
            else
                lowerB = gridOfCons(maxIx-1)
                upperB = gridOfCons(maxIx+1)
            end if
               
            maxDum = golden(lowerB, upperB, infoForMax, ForOptRoute, grids, valpol, params, tech, 1e-5_rk, argmaxDum)    

            ! Check that the new maximum is not much worse than the old one)
            ! Except whenn we are evaluating 'full' annuitization (and so full consumption of the remainder is optimal and the local searc will not reach the very top)
            if ((ixAnn .ne.grids%annPoints).and.(maxIx.ne.grids%cPoints)) then           
                if (maxDum.lt.candidateMax) then ! Try a greater tolderance
                    maxDum = golden(lowerB, upperB, infoForMax, ForOptRoute, grids, valpol, params, tech, 1e-20_rk, argmaxDum)    
                    if (maxDum.lt.candidateMax) then          
                        argmaxDum =candidateArgMax(1)
                        maxDum = candidateMax
                        if ((abs(candidateArgMax(1) - argmaxDum)).gt.10.0_rk) then ! If the difference in consumption space is less than 10.0 maybe don't worry about it
                            write(*,*) 'Local Supplement is worse than Candidate Max'
                            write(*,*) 'Annuity and optimal consumption indexes:', ixAnn, maxIx
                            write(*,*)  candidateArgMax(1), argmaxDum
                !            scratch = .true.
                        stop
                        end if                        
                    end if
                end if            
            end if
            
end subroutine

end subroutine getMax

!---------------------------------------------------------------------------------------------------------!
!---------------------------------------------------------------------------------------------------------!
! This is based on code given in Miranda and Fackler:
! "Applied Computaiontal Economics and Finance"

      real (kind=rk) function golden(a, b, infoForMax,  ForOptRoute, grids, valpol, params, tech, tol, argmax)
        implicit none

! Generic stuff
! Arguments
        real (kind=rk), intent (in) :: a
        real (kind=rk), intent (in) :: b
        type(typeGrids), intent(in) :: grids
        type(typeInfoForMax), intent(in) :: infoForMax
        type(typeForOptRout), intent(inout) :: ForOptRoute
        type(typeValPol), intent(in) :: valpol
        type(typeParams), intent(in) :: params
        type(typeTech), intent(in) :: tech    
        real (kind=rk) :: tol
        real (kind=rk), intent (out) :: argmax
   
! For program
        real (kind=rk) :: x1
        real (kind=rk) :: x2
        real (kind=rk) :: f1
        real (kind=rk) :: f2
        real (kind=rk) :: diff
        real (kind=rk) :: goldenalpha1
        real (kind=rk) :: goldenalpha2
        !real (kind=rk) :: tol
        parameter (goldenalpha2=0.618033988749895_rk)
        parameter (goldenalpha1=1.0_rk-goldenalpha2)
        !parameter (tol=1.d-30)

        x1 = a + goldenalpha1*(b-a)
        x2 = a + goldenalpha2*(b-a)

        ForOptRoute%leftOverA =  (ForOptRoute%AssetForLiqOrCon - x1)*(1+tech%r)      
        call getUtilPlusVal(infoForMax,  ForOptRoute, grids, valpol, params, tech, x1,f1)        

        ForOptRoute%leftOverA =  (ForOptRoute%AssetForLiqOrCon  -x2)*(1+tech%r)      
        call getUtilPlusVal(infoForMax, ForOptRoute, grids, valpol, params, tech, x2,f2)        

        diff = goldenalpha1*goldenalpha2*(b-a)
        do while (diff.gt.tol)
          diff = diff*goldenalpha2
          if (f2.lt.f1) then
            x2 = x1
            x1 = x1 - diff

            f2 = f1
            ForOptRoute%leftOverA =  (ForOptRoute%AssetForLiqOrCon  -x1)*(1+tech%r)      
            call getUtilPlusVal(infoForMax, ForOptRoute, grids, valpol, params, tech, x1,f1)        
          else
            x1 = x2
            x2 = x2 + diff

            f1 = f2
            ForOptRoute%leftOverA =  (ForOptRoute%AssetForLiqOrCon -x2)*(1+tech%r)      
            call getUtilPlusVal(infoForMax, ForOptRoute, grids, valpol, params, tech, x2,f2)        
          end if
        end do

        if (f2.gt.f1) then
          argmax = x2
          golden = f2
        else
          argmax = x1
          golden = f1
        end if

      end function golden



subroutine getUtilPlusVal(infoForMax, ForOptRoute, grids, valpol, params, tech, cons, utilPlusVal)
    implicit none

        type(typeGrids), intent(in) :: grids
        type(typeInfoForMax), intent(in) :: infoForMax
        type(typeForOptRout), intent(in) :: ForOptRoute
        type(typeValPol), intent(in) :: valpol
        type(typeParams), intent(in) :: params
        type(typeTech), intent(in) :: tech
        real (kind=rk), intent(in) :: cons 
        real (kind=rk), intent(out) :: utilPlusVal
        
        !For program
        real (kind=rk) :: ev1
        real (kind=rk) :: scratch 


        call getEV1_dev(infoForMax, ForOptRoute, grids, valpol%val, params, tech, ev1)            
        
        if (infoForMax%today.eq.1) then
            utilPlusVal = ev1 
        else if ((tech%ageAtSurvey + infoForMax%today -1).le.110) then           
            utilPlusVal = felicity(cons, params) + params%beta * ev1 
        else
            utilPlusVal= felicity(cons, params) 
        end if
        
                

end subroutine getUtilPlusVal

subroutine getEV1_dev(infoForMax,  ForOptRoute, grids, val, params, tech,  ev1, V1components)
    implicit none

        type(typeGrids), intent(in) :: grids
        type(typeInfoForMax), intent(in) :: infoForMax
        type(typeForOptRout), intent(in) :: ForOptRoute
        type(typeVal), intent(in) :: val
        type(typeParams), intent(in) :: params
        type(typeTech), intent(in) :: tech
        real (kind=rk), intent(out) :: ev1
        real (kind=rk), intent(out), optional :: V1components(2)
        
        ! For program
        real (kind=rk) :: leftOverA, v1 
        real (kind=rk) :: bothAssets

        

                    if (infoForMax%today .eq. (tech%retireAge - 1)) then ! In the last year of work (the year in which one can annuitize) income goes to whatever the new annuity stream will be
                            bothAssets = ForOptRoute%leftoverA
                            call linearinterp1(grids%A(:, infoForMax%tomorrow), val%vOldNoAccess(:, 1, ForOptRoute%ixAnn, infoForMax%tomorrow), grids%maxA, bothAssets, v1, 1, 1)
                    else     ! After retirement age one transits simply to the same income node as there are no earnings but the stream of annuity income is deterministic
                        if (infoForMax%solveT_SimF) then ! When we are solving, income is on the grid. When we are simulating, income tomorrow comes from infoForMax (fixed in maximiziation) rather than ForOptRoute (free in optimization)
                                  bothAssets = ForOptRoute%leftoverA
                                call linearinterp1(grids%A(:, infoForMax%tomorrow), val%vOldNoAccess(:, 1, infoForMax%ixY, infoForMax%tomorrow), grids%maxA, bothAssets, v1, 1, 1)
                             
                        else  ! When we simulating, income is not on the grid
                                bothAssets = ForOptRoute%leftoverA
                                    call linearinterp2(grids%A(:, infoForMax%tomorrow), grids%YOld(:, infoForMax%tomorrow), grids%maxA, grids%maxYO, bothAssets,  infoForMax%Y1, v1, val%vOldNoAccess(:,1, :, infoForMax%tomorrow))                                
                                
                        end if
                    end if   
                    !! Get expected value                    
                    ev1 = tech%surv(infoForMax%today) * v1 +  (1-tech%surv(infoForMax%today)) * 0.0_rk ! zero is the value being dead
                    
                 
                    ! Now do components for VSL if we're calling for it
                    if (present(V1components)) then    
                        V1components(1) = v1 ! if certain to live
                        V1components(1) = 0.0_rk ! if certain to die
                    end if
                
   
      
                




end subroutine getEV1_dev



!------------------------------------------------------------------------------------------------------------------------------------------------------------------!
!------------------------------------------------------------------------------------------------------------------------------------------------------------------!
      

    
    subroutine sumstats(series, years, numpeople, mean, sdev,  weights, logvals, topcode, quantiles)
        implicit none

!Arguments
        integer, intent (in) :: years
        integer, intent (in) :: numpeople
        real (kind=rk), intent (in) :: series(years, numpeople)
        
        real (kind=highprecintkind), intent (out), optional :: mean(years)
        real (kind=highprecintkind), intent (out), optional :: sdev(years)
        real (kind=highprecintkind), intent (out), optional :: quantiles(2, years)
        real (kind=rk), intent (in), optional :: weights(years, numpeople)
        logical, intent (in), optional :: logvals
        logical, intent (in), optional :: topcode
        
        
        
!Indexes
        integer :: ixyear, ixperson

!FOr programme
        integer :: haveweights
        integer :: ifail
        real (kind=rk) :: sumweights(years)
        real (kind=highprecintkind) :: xs_series(numpeople), quants(2)
        real (kind=highprecintkind) :: xs_weights(numpeople)
        real (kind=rk) :: seriesforprog(years, numpeople)

!Given that I am allowing a lot of optional arguments in the case of arguments that are not optional for the
!NAG routine I have to create duplicates that exist for the nag routine
        real (kind=rk) :: weightsfornag(years, numpeople)
        real (kind=highprecintkind) :: meanfornag(years)
        real (kind=highprecintkind) :: sdevfornag(years)
        real (kind=highprecintkind) :: skewfornag(years)
        real (kind=highprecintkind) :: kurtfornag(years)
        real (kind=highprecintkind) :: minfornag(years)
        real (kind=highprecintkind) :: maxfornag(years)
	    logical :: dologs, doTopCoding
	    real  (kind = rk) :: forTopcoding
	   	    real (kind=highprecintkind) :: rcomm(20)
            integer :: pn
            
        haveweights = 1
        weightsfornag = 1
        
        if (present(topcode)) then
            doTopCoding = topcode 
        else
            doTopCoding = .false.        
        end if

        if (present(weights)) then
          weightsfornag = weights
        end if

        if (present(logvals)) then !These have to be split into two conditions as logvals doesn't exist 
         if (logvals) then 
         dologs = .true.
         else
         dologs = .false.
         end if
         else
         dologs = .false.
        end if
         



        if (dologs) then
!The way I deal with zero or negative values is that I take the summary statistics ignoring these values
!Before taking the log, I replace these values with a positive number, to avoid an error
!But I set weights to be zero for these observations

          do ixperson = 1, numpeople, 1
            do ixyear = 1, years, 1

              if (series(ixyear,ixperson).le.0) then
                seriesforprog(ixyear, ixperson) = 0.0_rk
                weightsfornag(ixyear, ixperson) = 0.0_rk
              else
                seriesforprog(ixyear, ixperson) = log(series(ixyear,ixperson))
              end if

            end do
          end do
        else
          seriesforprog = series
        end if

!Want to get cross-sectional summary statistics
        do ixyear = 1, years, 1
          xs_series = real(seriesforprog(ixyear, :), highprecintkind)
          xs_weights = real(weightsfornag(ixyear, :), highprecintkind)
          
            ! If we're topcoding then topcode
            if (doTopCoding) then
                write(*,*) 'might need to check sumstats works with  new rks'
                stop
                ifail = 0
                call G01AMF(numpeople, xs_series, 1, 0.95_rk, forTopcoding, ifail)
                do ixperson = 1, numpeople, 1
                    if (xs_series(ixperson).gt.forTopcoding) xs_series(ixperson) = forTopcoding
                end do ! ixperson
            end if !   if (doTopCoding) then
                        
          if (sum(xs_weights).le.0) then                        
            
            meanfornag(ixyear) = -99999
            sdevfornag(ixyear) = -99999
            skewfornag(ixyear) = -99999
            kurtfornag(ixyear) = -99999
            minfornag(ixyear) = -99999
            maxfornag(ixyear) = -99999
          else
            haveweights = 1
            ifail = 1
            pn = 0
  
                if (numpeople.eq.1) then
                        meanfornag(ixyear) = xs_series(1)
                else
                            
                             call g01atf(numpeople, xs_series, haveweights,xs_weights, pn,meanfornag(ixyear), sdevfornag(ixyear), &
                             skewfornag(ixyear), kurtfornag(ixyear), minfornag(ixyear), maxfornag(ixyear),rcomm,ifail)

 
            if ((ifail.ne.2).and.(ifail.ne.0).and.(ifail.ne.71).and.(ifail.ne.72)) then
               write(*,*) 'Error in sumstats - NAG routine reporting a failure'            
            else if (ifail.eq.2) then !Only one valid case
             sdevfornag(ixyear)= -99999
             skewfornag(ixyear)= -99999
             kurtfornag(ixyear) = -99999
            end if 
            end if !if (numpeople.eq.1) then            
          end if
        end do

        if (present(mean)) mean = meanfornag
        if (present(sdev)) sdev = sdevfornag
if (present(quantiles)) then
    ! Now get quantiles
        quants(1) = 0.5
        quants(2) = 0.8
    do ixyear = 1, years, 1
         xs_series = real(seriesforprog(ixyear, :), highprecintkind)
         xs_weights = real(weightsfornag(ixyear, :), highprecintkind)
        call G01AMF(numpeople, xs_series, 2, quants, quantiles(:, ixyear), ifail)
    end do
end if

    end subroutine sumstats
!-----------------------------------------------------------------------------------------------------------!
!-----------------------------------------------------------------------------------------------------------!


    
    end module routines_specific